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ABSTRACT 

Despite the fact that the physics of the cosmic microwave background anisotropies is 
most naturally expressed in Fourier space, pixelised maps are almost always used in the 
analysis and simulation of microwave data. A complementary approach is investigated 
here, in which maps are used only in the visualisation of the data, and the tempera¬ 
ture anisotropies and polarization are only ever expressed in terms of their spherical 
multipoles. This approach has a number of advantages: there is no information loss 
(assuming a band-limited observation); deconvolution of asymmetric beam profiles 
and the temporal response of the instrument are naturally included; correlated noise 
can easily be taken into account, removing the need for additional ‘destriping’; polar¬ 
ization is also analysed in the same framework; and reliable estimates of the spherical 
multipoles of the sky and their errors are obtained directly for subsequent component 
separation and power spectrum estimation. The formalism required to analyse experi¬ 
ments which survey the full sky by scanning on circles is derived here, with particular 
emphasis on the Planck mission. A number of analytical results are obtained in the 
limit of simple scanning strategies. Although there are non-trivial computational ob¬ 
stacles to be overcome before the techniques described here can be implemented at 
high resolution, if these can be overcome the method should allow for a more robust 
return from the next generation of full-sky microwave background experiments. 

Key words: cosmic microwave background methods: analytical: - methods: nu¬ 
merical. 


1 INTRODUCTION 

The cosmic microwave background (CMB) represents the 
single most powerful probe of the early universe available 
to modern astronomy. The pioneering results of the Cosmic 
Background Explorer ( COBE ; Smoot et al. 1992) are be¬ 
ing extended by the current generation of ground-based and 
balloon-borne experiments (e.g. Scott et al. 1996; Tanaka et 
al. 1996; Netterfield et al. 1997; de Oliveira-Costa et al. 1998; 
Coble et al. 1999; de Bernardis et al. 2000; Hanany et al. 
2000; Wilson et al. 2000; Padin et al. 2001), but it is the up¬ 
coming satellite experiments that are set to revolutionise the 
field. By the end of 2002 the Microwave Anisotropy Probe ^ 
(MAP) will have observed the entire sky with a resolution 
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of ^ 0?21 and Planck ^, scheduled for launch in 2007, should 
have the capacity to produce high sensitivity, all-sky maps 
to a resolution of ~ 0? 1 in four of its ten frequency channels. 
However, it is not sufficient merely to obtain such extraordi¬ 
nary measurements; careful analysis of the data-sets is crit¬ 
ical if these missions are to fulfill their promised scientific 
goals. 

Most present techniques for CMB data analysis employ 
pixelised maps of the sky. Both Bond et al. (1999) and Bor- 
rill (1999) describe methods to create maximum-likelihood 
maps from time-ordered data and telescope pointing infor¬ 
mation, and these techniques have been successfully applied 
in the analysis of the BOOMERanG (de Bernardis et al. 
2000; Netterfield et al. 2001) and MAXIMA (Hanany et al. 
2000; Lee et al. 2001) observations. From the maps and their 
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uncertainties, the power spectrum of the microwave sky can 
be estimated in a number of different ways (Bond, Jaffe 
& Knox 1998; Borrill 1999; Oh, Spergel & Hinshaw 1999; 
Wandelt, Hivon & Gorski 2001), again employing maximum- 
likelihood techniques. These methods have a number of use¬ 
ful features: a map represents vast data compression rela¬ 
tive to the time-ordered data (e.g. Borrill 1999), but is still 
a sufficient statistic for cosmological parameter estimation; 
there are usually negligible pixel-pixel noise correlations in 
the beam-smoothed map; maps provide important ‘reality 
checks’, and can be inspected by eye; and unobserved or con¬ 
taminated regions of the sky can be simply removed from 
subsequent analysis by excluding the associated pixels (e.g. 
Oh et al. 1999). From a more practical point of view, the 
conventional map-making method is ‘tried and tested’ now. 
However, the use of pixelised maps during the data analy¬ 
sis also involves some compromises: the real microwave sky 
does not consist of a number of regions of uniform tem¬ 
perature (the prior hypothesis used to make a map), and 
so information is lost in map-making; subsequent steps in 
the analysis pipeline, such as component separation, cannot 
easily be performed efficiently in real space (Hobson et al. 
1998); there is no single optimal or obvious choice of pixeli- 
sation scheme; and deconvolution of the (asymmetric) beam 
profile and temporal response of instrument, and removal 
of scan-synchronous instrument effects require awkward ad¬ 
ditional processing during map-making (Delabrouillc 1998; 
Revenu et al. 2000; Wu et al. 2001). Another potential prob¬ 
lem that has received only limited attention thus-far (e.g. 
Tegmark & de Oliveira-Costa 2001) is how to create maps 
from polarization data (i.e. a map for each of the Stokes pa¬ 
rameters, or their gradient and curl components) - there is 
currently no robust algorithm for the treatment of several 
polarization-specific systematic effects that plague CMB po- 
larimetry experiments. 

It is with these points in mind that a complementary, 
harmonic method of data analysis is presented here. For full- 
sky surveys the data analysis process would proceed from the 
time-ordered data directly to the set of spherical multipole 
moments, which would take the place of a real-space map. 
For band-limited observations (effectively ensured by the fi¬ 
nite experimental beam), this harmonic reconstruction of 
the sky involves essentially no loss of information. The sep¬ 
aration of the various astrophysical components of the mi¬ 
crowave sky could then be performed quite naturally in mul- 
tipole space (e.g. Hobson et al. 1998), and, if necessary, the 
Galactic plane could be removed while retaining orthonor¬ 
mality of the basis set (Gorski 1994; Mortlock, Challinor 
& Hobson 2001; see also Section 5.2). Power spectrum es¬ 
timation could then proceed naturally from this point; for 
instance the efficient method presented by Oh et al. (1999) 
could be used, although some generalisation would be re¬ 
quired to remove the dependence on forward and inverse 
fast spherical transforms, which Oh et al. use to apply the 
inverse noise matrix efficiently in the space of multipoles and 
to reduce memory requirements. The multipole moments of 
the polarized components could be obtained using an almost 
identical formalism; obviously it would incur an additional 
overhead for each component to be reconstructed, but no 
further conceptual development would be required. 

We have endeavoured to present the harmonic method 
as an end-to-end solution for the processing of time-ordered 


data into power spectra. For such an analysis pipeline the 
benefits of clear error propagation from timelines to power 
spectra afforded by our harmonic methods are maximised. 
However, we do not pretend that all aspects of CMB analy¬ 
sis are best performed in the harmonic domain; instead we 
view harmonic methods as being complementary to standard 
map-based techniques. An obvious example where map- 
based processing is clearly desirable is the analysis of lo¬ 
calised foregrounds (or CMB features). In addition, the in¬ 
version from observations over only a fraction of the sky to 
the spherical multipoles rapidly becomes singular as the res¬ 
olution of the observations is increased. Although it may be 
possible to regularise the inversion in such cases (e.g. with 
prior information on the power spectrum), or to adopt a 
basis set adapted to the observed patch, the processing we 
describe in this paper is best considered in the context of 
full-sky observations. Even then, some aspects of the pro¬ 
cessing, such as cutting out contaminated regions close to 
the galactic plane for power spectrum estimation, require 
rather cumbersome operations if we insist on working di¬ 
rectly with the full sky spherical multipoles rather than a 
map synthesised from them. 

The principles of harmonic data analysis as described 
above are quite general, and could be applied to any full- 
sky CMB observations. However, it is particularly suited to 
experiments with a circular scanning strategy, such as the 
Planck mission. The data are obtained by a combination of 
rotations (i.e. the motion of the detector across the sky) and 
convolutions (of the beam with the true microwave sky, and 
this subsequent signal with the temporal response of the in¬ 
strument), and these operations are most simply expressed 
in the harmonic basis (Delabrouillc, Gorski & Hivon 1998a; 
Challinor et al. 2000; Wandelt & Gorski 2001). Hence we 
propose that the time-ordered data be transformed to the 
Fourier basis at the earliest opportunity: van Leeuwen et 
al. (2001) describe how to construct point source catalogues 
and calibrated ring-sets simultaneously, in which the data 
around each ring is represented in Fourier space. Some as¬ 
pects of this process are necessarily mission-specific, and, for 
the sake of generality, it is assumed here that it is possible 
to obtain calibrated ring-sets (and their errors) in this form. 
The focus of this paper is on the subsequent ‘map’-making: 
combining the ring-sets to obtain an estimate of the spheri¬ 
cal multipoles for the temperature anisotropies and polariza¬ 
tion. In this respect, we build on earlier work by Delabrouille 
et al. (1998a), who analysed the statistics of the power spec¬ 
trum of the Fourier modes on a single ring, and the relation 
with the power spectrum of the underlying temperature field 
on the sphere. More recently, Wandelt & Hansen (2001) have 
proposed the ‘ring-torus’ method for estimating the tem¬ 
perature power spectrum from the two-dimensional Fourier 
transform of data obtained on a set of rings with centres 
equally spaced on a small circle. While the general philoso¬ 
phy of the ‘ring-torus’ method coincides with that advocated 
here, the methodology presented here is somewhat more gen¬ 
eral. By only making use of one-dimensional Fourier data, 
we are able to deal with arbitrary ring-sets and to include 
instrument effects such as the inevitable variations in the 
scanning properties of the telescope (van Leeuwen et al. 
2001). Furthermore, since we do not proceed directly to the 
power spectrum, we can make full use of Fourier-based com¬ 
ponent separation algorithms (e.g. Hobson et al. 1998) to 
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deal carefully with foreground contamination. The cost of 
maintaining this generality is the increased computational 
requirements of our methods over those that exploit special 
symmetries of the survey. 

The general formalism is described in Section 2, which 
proceeds from a model of the instrument to the maximum- 
likelihood solution for the spherical multipoles of the tem¬ 
perature anisotropy and polarization, and their associated 
errors. The structure of the error covariance matrix, includ¬ 
ing polarization, is studied for some simple scan strategies in 
Section 3. Making use of some results given in Appendix B, 
we are able to make contact with several well-known analytic 
results for idealised experiments, previously obtained from 
arguments at the level of the map. In Section 4 we discuss 
a number of important issues that arise during the map¬ 
making phase, including correcting for scan-synchronous in¬ 
strument effects, the overall calibration of the instrument 
to external standards, and the treatment of low frequency 
noise. The latter discussion includes a novel method for deal¬ 
ing with uncertainties in the low frequency noise power spec¬ 
trum or insensitivity of the experiment to the monopole. 
Section 5 reviews the subsequent processing of the frequency 
maps, including component separation and power spectrum 
estimation, within the context of the harmonic data model. 
Finally, we conclude in Section 6 by reviewing the relative 
merits of the harmonic method, and suggest directions for 
future development. 


2 THE HARMONIC DATA MODEL 

The basis of the harmonic data model is the relationship be¬ 
tween the one-dimensional Fourier representation of the ring 
data and the spherical multipole coefficients of the underly¬ 
ing sky. In the absence of instrument noise and reconstruc¬ 
tion errors in the estimation of the Fourier ring data, the 
relation between the data and the sky is determined by the 
scanning strategy, the point-spread function or beam of the 
telescope, and the impulse response of the instrument in the 
time domain. In this section, we set up a detailed, but gen¬ 
eral, model of the instrument and scanning strategy which 
defines the contribution to the data coming from the sky. 
We then include instrument noise and reconstruction errors 
to give the maximum-likelihood solution for the spherical 
multipoles of the sky, and the covariance of their errors. 


2.1 The microwave sky 

We describe the microwave sky in terms of spherical multi¬ 
poles on some fixed basis. The brightness in total intensity 
/(e; v) from the sky when observed along direction e at fre¬ 
quency v is represented as 

I(e; v) = ^ a* lm )(v)Y ( im)(e), (1) 

lm 

where the sum is over integers 1^0, and \m\ ^ l. The (par¬ 
tial) polarization of the sky is described by Stokes (bright¬ 
ness) parameters Q{e;v), U(e;v), and V{e;v). We define 
the Q and U Stokes parameters on a basis which forms 
a right-handed triad with the incoming radiation direction 


—e: In a spherical polar coordinate system, we use the ba¬ 
sis {erg,— cr^} with oy = e. The V brightness is a scalar 
function, so can be expanded in multipoles as 

V(e-u) = 0(i m) (y)Y(i m ) (e). (2) 

lm 

It is convenient to combine the Q and U brightnesses into a 
symmetric, trace-free second-rank tensor: 

V ab {e;v) = i[Q(e; u)(a e (g> erg - 0 af) 

— U{e\v){<re ® erg)], (3) 

which can then be expanded in the transverse, trace-free ten¬ 
sor harmonics (Kamionkowski, Kosowsky & Stebbins 1997) 
as 

Table; v) = + 0 (r m )(^)F(^ m ) oi) (e)]. (4) 

lm 

Here, the sum is over l ^ 2, and \m\ ^ l. The super¬ 
scripts G (for gradient, often called electric) and C (for 
curl, often called magnetic) refer to the two types of trans¬ 
verse, trace-free harmonics. All multipoles satisfy offf n )(y) = 
(—1 ) m ffl(j_ m )(i / ) since the brightnesses are real and we have 
adopted the Condon-Shortley phase for the spherical har¬ 
monics. 

In this paper we assume that the af lm fv) contain all 
astrophysical components. In particular, we include the con¬ 
tribution of unresolved extra-Galactic radio sources (point 
sources). The point source catalogue is an important deliv¬ 
erable product of the Planck mission; an algorithm for si¬ 
multaneously constructing the bright point source catalogue 
and calibrating the instrument from data ordered by phase 
on rings is described by van Leeuwen et al. (2001). Here, we 
assume that the contribution of those bright point sources 
identifiable on individual rings has not been removed from 
the data, although we leave open the question of whether 
this assumption is optimal in the context of the harmonic 
data model. The steps described in this paper for recon¬ 
structing the multipoles of the sky from the Fourier modes 
of the ring data are largely independent of whether the 
data includes the bright point source contribution or not. 
Of course, if the bright point sources are removed from the 
data prior to reconstructing the sky, our solution excludes 
the contribution from those identified point sources. One 
practical problem with leaving bright point sources in the 
ring-sets is that the sky must be represented by a larger 
number of multipoles to avoid biasing the estimates of the 
lower multipoles. Furthermore, it is then essential that the 
point sources are removed during the component separation 
phase, before any cosmological analysis of the maps (e.g. 
Vielva et al. 2001). Removing the brightest point sources at 
the level of the ring-sets leads to a modest increase in the 
complexity of the pipeline leading from time-ordered data to 
calibrated ring-sets. In particular, although a given source 
will only give rise to a prominent feature on a small subset 
of rings, the source’s contribution must be removed consis¬ 
tently from every ring to ensure that the remaining data 
on each ring is derived from the same underlying sky. Point 
sources that are too faint to be identifiable on any single 
ring may still be detectable statistically when the rings are 
combined into the best-fitting set of multipoles (or map). 
Statistical detection of faint point sources is easily imple- 
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mented during the component separation phase (Hobson et 
al. 1999), and does not require the same pre-processing steps 
as removal of bright point sources (Vielva et al. 2001). 

2.2 Ring data 

A CMB mission with a circular scanning strategy can be de¬ 
scribed in terms of JVd sets of N r rings on the sky, each set 
referring to a given detector (see Fig. 1). For all detectors, 
each ring is covered N s times by spinning the instrument 
about the axis of the ring. In the case of the Planck mis¬ 
sion, the satellite will rotate N s ~ 60 times about it axis 
before moving on to a new pointing. The rth ring is speci¬ 
fied by its axis pointing, (0 r ,<p r ), which coincides with the 
average pointing of the spin axis during the N s revolutions. 
The axis pointing of a ring is the same for all detectors on 
the instrument. Note that the concept of a ring only makes 
sense if any precession of the spin axis (due to initial mis¬ 
alignment of the spin axis with the principal directions of 
the inertial tensor, or external torques) can be engineered 
to be insignificant compared to the smallest beam width. 
The opening angle of the rth ring for the dth detector is 
denoted by a r d- It is defined as the average (over the N s 
revolutions) angle between the spin axis and the nominal 
main beam direction on the sky of the appropriate horn in 
the focal plane. Not only does a r d depend on the detec¬ 
tor, but it also differs between rings for a given detector. 
The dependence on the detector is determined by the focal 
plane geometry. For low frequency, polarization-sensitive in¬ 
struments (e.g. HEMT receivers) there will typically be two 
detectors (measuring nearly orthogonal polarization states) 
on a given horn. In such an arrangement, the ring opening 
angle will be the same for the two detectors, but will differ 
between horns. Similar comments apply to the current state 
of the art designs for high frequency bolometer instruments, 
where two orthogonal polarization sensitive bolometers are 
placed on a single horn. The dependence of a r d within a 
given detector’s ring-set arises from slow drifts of the spin 
axis relative to the instrument optics as consumables are de¬ 
pleted during the mission (thus changing the inertia tensor 
of the instrument). 

The position around a given ring is specified by ip, and 
is measured so that ip = 0 is the most southerly intersection 
of the ring and the great circle contaning the 2 -direction and 
the spin axis [see Fig. 1; for 6 r = 0 (7r) we take ip = 0 to 
lie in the x-z plane with positive (negative) x\. The times at 
which different detectors pass through ip = 0 is dependent 
on the focal plane geometry. For each horn in the focal plane 
we define a constant reference configuration where the nom¬ 
inal main beam direction is aligned with the <j~ direction of 
the Cartesian frame used to define the spherical multipoles 
for the sky, and the 2 -axis of the instrument reference sys¬ 
tem (van Leeuwen et al. 2001) lies in the plane normal to cr y , 
with a negative projection onto cr x . (The instrument refer¬ 
ence system, which is fixed relative to the instrument optics, 
can be chosen so that its 2 -axis almost coincides with the 
nutation-averaged spin axis. Variations in the inertia ten¬ 
sor during the mission prevent a constant alignment of the 
spin and 2 -axes.) The instrument configuration when the 
dth detector takes data at angle ip on the rth ring is ob¬ 
tained by rotating the instrument from the appropriate horn 
reference configuration. The appropriate rotation is given 


by the composition-! D{cp r , 6 r , ip)D(0, a r d, n r d). The rota¬ 
tion D(0, a r d, K r d) accounts for focal plane rotation (which 
arises from misalignment of the spin axis and the 2 -axis of 
the instrument) by: (i) first rotating the spin axis into the x- 
2 plane of the sky coordinate system while leaving the main 
beam direction along the 2 -direction on the sky; and (ii) tak¬ 
ing the spin axis onto the 2 -direction on the sky by rotating 
in the plane contaning the spin axis and the main beam di¬ 
rection. The angle K r d thus measures the angle between two 
planes, both containing the main beam direction, with one 
including the spin axis and the other including the 2 -axis of 
the instrument reference system. The rotation D{(p r ,6 r ,ip) 
takes the spin axis onto the axis of the rth ring, and the 
main beam to the position ip around the ring. Note that the 
horn reference configuration is defined by the instrument ref¬ 
erence system rather then the spin axis. This ensures that 
the beam patterns in the given horn reference configuration 
are independent of variations in the inertia tensor. 

The contribution of the (time-independent) sky signal 
to the time-ordered data from a given detector will be peri¬ 
odic in ip, so the data can be co-added in bins of ip with es¬ 
sentially no loss of useful information. Producing this phase- 
ordered data from the time streams is a non-trivial task, 
since typically the scan velocity drifts during the N s scans 
of the circle: van Leeuwen et al. (2001) detail one possible 
scheme for reconstructing the phase-ordered data and asso¬ 
ciated errors. From the phase-ordered data t r d{ip), one can 
estimate the Fourier coefficients in the expansion 

trd{ip) = rdn )e n ^, (5) 

n 

and the covariance of their errors N^^^r'd'n') = 

(At( r d n )At( r , d , n ,\). Here, angle brackets denote the expec¬ 
tation value and At( r dn) — f(r-dn) (f(rdn))- The t( r dn) aie 
the primary data objects in the harmonic reconstruction of 
the microwave sky. 

2.3 The instrument response 

We assume that the detector responses are linear function¬ 
als of the sky in the absence of instrument noise. At any 
instant, the power seen by a given detector (due to the sky) 
is given by convolving the sky with the appropriate beam 
and spectral transmission. A single time-ordered data point 
is then obtained by convolving this power with the detector’s 
temporal response. 

2.3.1 Beam patterns 

It is convenient to define the beam profile when the in¬ 
strument is in the appropriate horn reference configura¬ 
tion, so that the beams can be assumed not to vary 
through the mission. A method for describing arbitrary, 
polarized beam patterns in terms of multipoles b^ lrn ^ was 

§ Our convention for the Euler angles a, 0 and 7 are such that 
the rotation D(a,/ 3 , 7 ) actively rotates by 7 about rr z . followed 
by 0 about cr y , and finally by a about <r» again. All rotations are 
right-handed. See Brink & Satchler (1993), whose conventions we 
follow, for a discussion of the several alternatives that appear in 
the literature. 
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Figure 1. Geometry of the rth ring for the <ith detector. The 
average direction of the spin axis over the N s revolutions that 
make up the data on a single ring has polar angle 6 r and azimuth 
(f>r , and is the same for all detectors. The ring opening angle a r d 
depends on the detector and the ring (due to variation in the in¬ 
ertia tensor during the mission). For a given ring and detector, ^ 
measures the angle around the ring. The point -if) = 0 is defined as 
the most southerly intersection of the ring and the plane contain¬ 
ing the ring axis and the ^-direction. Note that ^ = 0 is attained 
at different times for different detectors - the phase offsets being 
determined by the focal plane geometry. 

introduced by Challinor et al. (2000). The beam pattern 
is first described in terms of effective Stokes parameters 
{Id(e;iy),Qd(e;jy),Ud(e;u),Vd(e]iy)}, which are defined on 
the same {cr#, — cr^} basis as the sky. These effective Stokes 
parameters can be expressed in terms of the far-field radi¬ 
ation pattern of the detector (Challinor et al. 2000). The 
Qd(e;z/) and C7d(e;z/) can be combined into a linear polar¬ 
ization tensor for the beam, #5 b (e;z/), as in equation (3). 
The resultant power per unit frequency interval, dWd/du , 
can then be written in basis-independent form as the inte¬ 
gral 

Wd = A eS ,d J[I(e-u)I d (e-u)-V{e;u)Vd{e-u) 

+ 2Vab{e-,v)B°i b {e-,v)}dn, (6) 

where A e ff,d is the effective area (assumed independent of 
frequency). The directivity I d {e\v) is a scalar field on the 
sphere, so can be expanded in terms of multipoles 6 d(im)(^)> 
as in equation (3). We impose the constraint 6 ^( 00 ) (t) = 
so that /d(e; u) integrates to unity over the sphere. In 
a similar manner, Vd(e; v) can be represented by multipoles 
Finally, Bd b (e-,u) can be expanded in the trans¬ 
verse, trace-free tensor harmonics with multipoles ( v ) 

and bf (lm) H)- 

To obtain the power per unit frequency interval inci¬ 
dent on the spectral filter chain for the dth detector, when 
the nominal main beam is at angle ■)/> on the rth ring, we 
must rotate the beam pattern by 9 r ,ip)D{ 0, aw, fiw) 


before performing the convolution in equation (6). This ro¬ 
tation is easily performed in multipole space to give the re¬ 
sult (Challinor et al. 2000) 

-r; W rd {i >) = Aeff.d ^2 [ a O*rn){v)bd(im’){v) 

Pimm'm" 

X-Dmm" l{i)Dm" m' (0? &rd, ^rd)]) (7) 

where, for later convenience, we have defined 



+ 25Gbd( lm fv)+28cbd( lm )(v), ( 8 ) 

with bf (lm) {is) = bf (lm) {v) = 0 for l < 2. The D l mrn , (a, /3, 7 ) 
appearing in equation (7) are Wigner’s D-matrices; our con¬ 
ventions follow Brink & Satchler (1993). The dependence of 
OL'(«,A 7 ) on a and 7 goes as 

D l mm , (a, / 3 , 7 ) = e- im “<4w (9) 

where the d l mrn , {(3) are the reduced rotation matrices. 

2.3.2 Spectral filtering 

The power available at the detection element (e.g. bolome¬ 
ters for high frequency instruments) is given by integrating 
dW r d(if)/di> against the spectral transmission Vd{v) of the 
filter chain. We shall assume that the filters on all detec¬ 
tors in a given frequency band are identical, and have trans¬ 
mission properties that are accurately known from calibra¬ 
tions on the ground. Map-making from multi-frequency data 
is usually performed on a band-by-band basis (although it 
is essential to use all the frequency bands to calibrate the 
pointing, e.g. van Leeuwen et al. 2001). Given that Fourier 
data from only detectors in the same frequency band are 
combined in such an analysis, we need only consider a single 
band, and henceforth assume that v d {v) has no dependence 
on the detector, i.e. v d {v) — v{v). 

On integrating equation (7) against v{v), the de¬ 
pendence on the sky enters through the integral 
/ a fim){ v )b d (irn'){ l ') v { l ')dv. For the subsequent analysis, it 
is important that this integral can be factored into a part 
describing the sky carrying only l and m indices, and a part 
describing the beam which has only l and m! indices. Here 
we assume that the approximate factorisation 

J a(lm)(v)bd(lm')(v)v(v)dl' « o) 

X J dv, (10) 

where uo is the central frequency in the band, is accurate to 
better than the noise level, so that we can describe the sky 
by the frequency integrated multipoles 

d(im) = J Him){v)v(y)dv. ( 11 ) 

The largest source of error resulting from this factorisation is 
expected to be due to the variation of the width of the main 
beam across the band. In Appendix A we compute the sys¬ 
tematic error on the recovered power spectrum of the CMB 
temperature anisotropies resulting from this variation, as¬ 
suming an axisymmetric, diffraction-limited Gaussian beam, 
uniform sky coverage, and no foreground contaminants. For 
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the Planck High Frequency Instrument (HFI), the bias may 
be non-negligible in the 100 GHz channel between l = 1000 
and l = 1500, in which case refinements to equation (10) 
will be required. Here, we adopt the factorisation in equa¬ 
tion (10), in which case the power available at the detection 
element, W rd (ip), can be written as 

W r d(lp) = A e ff } d ^ (<Pr, @r, Ip) 

Pimm' 

X Brd(lm')]- (12) 


For later convenience, we have introduced the notation 

^rd(lm') — ^ ' Dm'm" (Q> CXrd, ^'rd')bd(lm") (^o) (13) 

m" 

for the multipoles of the dth beam after rotation by 
D(0,a r d,K r d). Note that the r-dependence of de¬ 

rives only from the (slow) variation of opening angle a r d 
and focal plane rotation k r d through the mission. 

2.3.3 Temporal response 

The available power, W rd [ip{t')\, is convolved with the tem¬ 
poral response function of the detector, hd(t,t'), to gener¬ 
ate a single time-ordered datum at time t. We decompose 
hd(t,t') into a stationary part, hd(t — t 1 ), and a stochastic 
part, A hd(t,t'), which we assume to have stationary statis¬ 
tical properties. The stochastic part, which derives from e.g. 
random fluctuations in amplifier gains, will be considered as 
part of the instrument noise in Section 2.4. Here, we con¬ 
centrate on the stationary impulse response function h d {t). 
Typically, this will contain several instrument artifacts, such 
as the effect of finite sampling and non-zero detector time 
constants, which are convolved together to form the total 
impulse response. For example, if the detector time constant 
is Td, and samples are taken by integrating over time inter¬ 
vals A d and are assigned to the midpoint of the interval, 
then hd{t) = h Tdtd (t) * hA d ,d(t), where the detector impulse 
response 

h Td ,d(t) = G d Td~ 1 e~ t/rd Q(t/Td), (14) 

and the sampling impulse response 

hA d ,d(t) = A<z _1 0(l + 2t/A d )0(l - 2t/A d ). (15) 

Here, ©(*) is the Heaviside unit step function, and G d is the 
average gain over the mission, which includes both amplifier 
gain and the quantum efficiency of the detectors. 

One of the advantages of describing the data on rings 
in the Fourier domain is that the convolution with the tem¬ 
poral response of the instrument reduces to a simple prod¬ 
uct (Dclabrouille et al. 1998a). Denoting the contribution of 
the signal to the phase-ordered data by s rd {ip), we have 

/ OO 

W r d{tp + OJrt')hd(-t') dt', (16) 

-CO 

where u> r = dip/dt is the (average) spin rate on the rth ring. 
We have assumed that the support of h d (t) is sufficiently 
compact that the only contribution to s rd {ip) is from the 
available power W rd on the same ring. We have further as¬ 
sumed that any variation of the spin rate about the average 
is negligible for the purpose of performing the convolution, 


so that the convolution of the available power with the im¬ 
pulse response is periodic in ip. Inserting equation (12) into 
equation (16), and expanding D l mm , {<p r , 6 r , ip) as in equa¬ 
tion (9), we find 

Srd{lp) = A e ff ? d y ' [h(; m )D mm / (4>r, 9r, Ip) 

Pimm' 

X Hd(m'u>r)Br d (i m >)], (17) 


where Hd(co r ) is the Fourier transform of the impulse re¬ 
sponse: 


Hd{yOr) 



h d (t)e~ iuJrt dt. 


(18) 


For the impulse responses in equations (14) and (15), we 
have 


Hdip} r) 


G d sinc(aj r A d /2) 
1 + iuj r Td 


(19) 


Finally, we can extract from equation (17) the contribution 
of the signal to the nth Fourier coefficient of the phase- 
ordered data. Complex conjugating and using equation (9), 
we find 


S(rdn) = A eS }d afl m) dl nrl (e r )e rnrl,r H d (nUJ r )Brd(ln )], (20) 

Plm 


where the sum is over l ^ |n| and |m| ^ l. Note that the 
effect of the impulse response can also be described in terms 
of an effective (ring-dependent) beam with multipoles 

b rd(fm)( U ) = X! [ D l mm'(-Krd,-a r d,0)Hd(m'uj r ) 
m'm" 

X Dm'm" (0? CXrdi ^rd)^d{l m") Ml-(21) 


The effective beam will be generally not be axisymmetric 
as a consequence of the smearing along the scan direction 
induced by the impulse response of the instrument. The con¬ 
cept of an effective beam becomes particularly significant if 
the beam has to be calibrated from the observations of e.g. 
bright point source transits, as described by van Leeuwen et 
al. (2001). In such cases, it is the effective beam that can 
be reconstructed most directly. A potential difficulty that 
would need to be addressed in such an approach is the (slow) 
variation of effective beam with ring due to variations in the 
angles a rd and K rd . In practice, it may be simpler to decon¬ 
volve the temporal response of the instrument prior to con¬ 
struction of the phase-ordered data. In this case, H d (nu> r ) 
can be omitted from equation (20). 


2.4 Gaussian random noise and reconstruction 
errors 

The Fourier modes, t( rdn ), that we construct from the phase- 
ordered data differ from the S( rdn ) of equation (20) because 
of a number of sources of error. For example, stochastic- 
ity in the instrument temporal response, such as that due 
to thermal noise in the detectors and amplifiers, has not 
been included in equation (20). Neither have we included 
the effects of photon (shot) noise, or any detector offsets. 
Furthermore, we have replaced the spin axis pointing by 
its average over the N s spin periods, ignoring the (albeit 
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small) effects of the nutation of the instrument. These ran¬ 
dom errors tend to be suppressed in the phase-ordered data 
due to the averaging over N s periods, except for fluctua¬ 
tions at temporal frequencies synchronous with the spin fre¬ 
quency. Fluctuations below the spin frequency remain in the 
phase-ordered data predominantly as an offset on the ring. 
In the usual map-making paradigm (Bond et al. 1999; Borrill 

1999) , low frequency (1//) noise must be carefully accounted 
for; failure to do so leads to highly correlated errors in the 
reconstructed map which typically appear as stripes (e.g. 
Delabrouille 1998; Revenu et al. 2000). It is difficult to ac¬ 
commodate low frequency noise power in map-making from 
the time-ordered data since its inclusion spoils the sparse na¬ 
ture of the time-time noise correlation matrix. Approximate 
methods to circumvent this problem for experiments that 
scan on rings already exist (Delabrouille 1998; Revenu et al. 

2000) ; these ‘destriping’ methods will be discussed further 
in Section 4.1. 

There are also a number of potential systematic errors 
which contaminate the t( r dn), and, if left unaccounted for, 
will lead to bias in the reconstructed multipoles of the sky. 
Keeping the numerous systematics under control is essential 
to maintaining the integrity of the final data products. A sig¬ 
nificant merit of the harmonic data model is that the model 
of the instrument, equation (20), is sufficiently complete to 
include a number of systematic effects which are not so nat¬ 
urally included in the standard map-making techniques. For 
example, the harmonic model can seamlessly accommodate 
any prior knowledge we have of asymmetries in the beam 
profiles, cross-polar contamination for the polarized detec¬ 
tors, sidelobe features leading to ‘straylight’ entering the 
instrument, and detector time constants and data sampling 
rate. For some systematics there will be incomplete (or even 
no) prior knowledge of essential parameters, in which case 
these parameters must be estimated from the data them¬ 
selves during the reconstruction of the sky. The harmonic 
model also appears well-suited to handling such systematics 
iteratively; see Section 4.2. 

To include the effect of Gaussian random noise with a 
known (or estimated) power spectrum in the reconstruction 
of the sky from the Fourier data on rings, t( r dn ), we require 
the covariance matrix N^d^^r'd'n')- This matrix was calcu¬ 
lated by Delabrouille at al. (1998a) for a single detector and 
ring (d = d! and r = r'). Here we extend their result to in¬ 
clude the effects of a finite sampling rate (see also Janssen et 
al. 1996), and to allow for r ^ r'. Assuming instrument noise 
dominates the random errors Att r dn)j there will be negligi¬ 
ble correlation between the errors on different detectors, so 
we need only consider d = d!. 

We assume that the noise contribution, n<j(t), to the 
time stream of the dth detector is a real, stationary random 
process with zero mean, and power spectrum Md[co). The 
power spectrum is the Fourier transform of the noise auto¬ 
correlation function Cd(r) = {nd(t + r)nd(t)): 

/ OO 

C , d(r)e _I " T dr. (22) 

-OO 

The propagation of the noise rid(t) to the errors A t( rdn ) 
depends on the exact procedure for estimating the Fourier 
modes t( r d n ) from the time-ordered data (van Leeuwen et al. 

2001). However, to a first approximation, the errors At( rdn ) 
are obtained by convolving rid(t) with the impulse response 


of the sampler, hA d ,d(t), mapping the portion of this func¬ 
tion covering the observation interval for the rth ring onto 
the angle i/j around the ring, and finally extracting the 
Fourier coefficients of the resultant signal: 

r t+ 

LO I 

A t.(rdn) = 2 ^ ft-A d,d * n d (f) exp[-muv(f - t~d )] d t. (23) 

Here, t~ d is the time when the dth detector is first pointing 
to ip — 0 on the rth ring, tf d is the end of the observation 
interval for that ring: tf d = t~ d + 2nN s /u} r , and we have 
ignored any variation of the spin rate during the N s revo¬ 
lutions. It is now straightforward to write the covariance of 
the errors in terms of the noise power spectrum. For the case 
of scanning missions such as Planck, the number of revolu¬ 
tions N s 1 is sufficiently large that negligible correlations 
remain between different Fourier modes. In this case we find 

1 f 00 

N( r dn)(r'd'n') ^dd'^nn'’^ / {A/*d(^)sinC (caAd/2) 

x sinc 2 [(u;/t<; r — n)7rAa] 
x exp[»27rAf s (r — r')ui/ ca r ]} du>, (24) 

where we have further ignored any variation of the 
average spin rate, uv, between rings. Provided that 
Nd{a>)smc 2 {uiAd/2) varies slowly compared to the rest of 
the integrand in equation (24) around u> = nw r , we can re¬ 
place it by its value there. The remaining integral vanishes 
unless r = r' , so that for slowly varying power spectra the 
covariance matrix is fully diagonal: 

N^ r dn)(r'd'n') — S rr ' $dd' ^nn' 

x -^-A/d(nca 7 -)sinc 2 (nw r Ad/2). (25) 

27TiV s 

In the presence of a significant low frequency noise com¬ 
ponent (such as 1 // noise in the electronics) joining onto 
white noise at a knee frequency Wknee (e.g. Delabrouille 
1998), the assumption of a slowly varying noise power spec¬ 
trum may not hold for small |n|. Typically the spin rate will 
be chosen so that Wknee u> r , in which case equation (25) 
will be valid for n ^ 0. However, for n = 0 we can ex¬ 
pect correlations between rings unless ui r /N s <C u>knee -C ui r - 
(For the Planck HFI, the nominal Wknee < 0.06rads -1 , and 
ui r « 0.10rads -1 , so these conditions are likely to be sat¬ 
isfied.) Even in the presence of such correlations, the noise 
covariance matrix is still block diagonal in the harmonic rep¬ 
resentation, in contrast to the time-time covariance matrix 
which has significant off-diagonal terms due to the extended 
support of the noise auto-correlation function, Cd(r). 

2.5 Maximum-likelihood solution 

Our model for the ring Fourier data t( r d n ) is now t( r dn) = 
S(rdn) + At( r dn), where the signal S( rd n) is modelled by equa¬ 
tion (20). It is convenient to write this relation in the form 

t'(rdn) E -A(rdn)(Plm)Q J (lm) “1“ ^^(rdn) ? (26) 

Plm 

where the coupling matrix is 

A( r dn)(Plm) = ^4eff ,^< 4^71 (fir ) el "’ ^ Hd (nUr )Bf d (ln) ■ (27) 

We define A( rd n)(pim) to be zero for l < |n|, since the t( rdn ) 
only depend on the df lrn ^ with l ^ |n|. This structure can 
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be exploited to derive an unbiased (although not minimal 
variance) estimator for the hp m ), which can be computed 
in 0(lm ax) operations. Here, Z ma x is the maximum Z that we 
retain during the inversion. (The experimental beam limits 
imax even if the underlying sky is not band-limited.) The 
estimate is derived by working down from the maximum 
Fourier mode, n max , and performing a regularised inversion 
of a subset of the data to solve for those on which the 
data subset depends but which have not been determined at 
a previous step. Full details of this estimator will be given 
in a future paper (Mortlock et al. in preparation). 

In this paper we shall just give the formal maximum- 
likelihood solution to equation (26). We assume that we 
are attempting to solve only for the which form the 

components of a vector a. In practice, it may also be de¬ 
sirable to attempt to solve for several instrument parame¬ 
ters which are unknown from earlier calibrations, but which 
represent significant systematic effects. Such parameters are 
easily included in the maximum-likelihood formalism, but 
unless they influence the data linearly their inclusion would 
prevent us from being able to locate the maximum-likelihood 
multipoles analytically. Given the size of the problem, hav¬ 
ing to locate the maximum of the likelihood by a numerical 
search in model space is clearly undesirable. The treatment 
of systematics with unknown parameters is probably best 
performed iteratively (e.g. Delabrouille, Gispert & Puget 
1998b; see also Section 4.2). 

Assuming Gaussian noise, the maximum-likelihood es¬ 
timate, a, for the true sky, a, is 

o-,== (A t N _1 A) _1 A t N _1 t, (28) 

where t is the vector of Fourier modes t( r( j n ), A is the ma¬ 
trix of coupling coefficients A^ r d n )i r 'd'n')i and N is the noise 
covariance matrix with components N( r dn)(r'd'n')- The ^ op¬ 
eration in equation (28) denotes Hermitian conjugation. The 
maximum-likelihood solution is the optimal, unbiased esti¬ 
mate of the sky. Its covariance matrix is 

C = ((a — a)(a — a) f ) = (A t N _1 A) _1 , (29) 

where the average is over noise realisations. A brute force 
evaluation of a requires 0(Zm ax ) operations and 0(Zm ax ) 
storage which is impossible at the moment for high resolu¬ 
tion experiments. Conjugate gradient techniques (e.g. Press 
et al. 1992) would reduce the operations count to 0(Nil^ ax ), 
where N-, is the number of iterations required, by removing 
the matrix multiplies and inversion. Note that the block- 
diagonal structure of N (equation 24) reduces the cost of 
computing N _1 to 0(Zm ax ) operations, or fewer if correla¬ 
tions between rings are confined to the n = 0 modes. Keep¬ 
ing Ni small requires a careful choice of preconditioner, i.e. 
a Hermitian, positive-definite matrix which is both a good 
approximation to the inverse covariance matrix C _1 and is 
easy to invert. We demonstrate in Section 3 that for simple 
scan strategies such matrices can easily be found (see also 
Oh et al. 1999). For the case where all rings are at similar lat¬ 
itude, an approximation to C _1 can be computed in 0(Zm ax ) 
operations and requires only 0(Z alax ) storage. This precon¬ 
ditioner is block-diagonal and can be inverted in 0(Z((, ax ) 
operations. Note, however, that the conjugate gradient eval¬ 
uation of equation (29) still requires 0(Z^ ax ) storage for the 
non-sparse matrix A. The application of the conjugate gra¬ 
dient method, as well other iterative techniques, to equa¬ 


tion (29) will be explored numerically in Mortlock et al. (in 
preparation). 

In writing down the maximum-likelihood solution we 
have assumed that A^N *A is invertible, so that the solu¬ 
tion is unique. For some scan strategies and instrument ge¬ 
ometries, A^N X A may become numerically singular due to 
incomplete coverage of the sky. With incomplete coverage, 
and Z max sufficiently large, there may exist sets of multipoles 
which produce a temperature or polarization field which is 
localised to machine precision in the regions of the sky that 
are not covered (Mortlock et al. 2001). In such cases, solv¬ 
ing for the multipoles must proceed via singular value tech¬ 
niques, and we obtain no constraint on the part of a which 
lies in the null space of the coupling matrix A. If map-making 
were performed in pixel space, similar problems would arise 
when attempting to estimate the multipoles from the pix- 
elised map. As an alternative to singular-value techniques, 
we could further regularise the inversion by Wiener filter¬ 
ing. With a Gaussian prior for the signal a, with covariance 
S = (aa 1 ), the maximum of the posterior probability gives 
the Wiener-filtered reconstruction 

a = (S' 1 + A t N~ 1 A) _1 A t N _1 t. (30) 

In this manner, the inversion of the singular matrix A 1 N -1 A 
is regularised by the inverse signal covariance matrix S -1 . In 
practice, the microwave sky is only poorly approximated as 
a Gaussian random process, particularly in those frequency 
channels where the CMB is not dominant. In addition, we 
will only have limited knowledge of the signal covariance 
matrix S either from other observations, or a preliminary, 
approximate analysis of the data. However, Wiener filtering 
has proved to be a useful technique in component separation 
where similar objections to the use of Gaussian priors hold 
(see e.g. Hobson et al. 1998 and references therein), and 
the same should be true for harmonic ‘map’-making. Note 
that Wiener filtering produces a biased estimate of the sky; 
other linear filters have been proposed to cirumvent this 
problem, but generally produce noisier, though unbiased, 
maps (e.g. Tegmark & Efstathiou 1996). Equation (30) can 
also be solved with conjugate gradient techniques; in this 
case the preconditioner should include some simplified form 
of the signal covariance in addition to an approximation to 
A + N^A . 

In our discussion so far we have assumed that full beam 
deconvolution is performed during the ‘map’-making stage. 
In practice, deconvolving the beam completely may be un¬ 
desirable for several reasons. Firstly, the matrix A 1 N _1 A is 
likely to be more ill-conditioned since the errors on the re¬ 
constructed multipoles must grow large as Z approaches the 
(inverse) beam scale. (In a conjugate gradients solution, the 
preconditioner can partly take this effect into account.) Sec¬ 
ondly, as discussed in Section 2.3.2 and Appendix A, varia¬ 
tion of the width of the main beam across the spectral band 
of a channel is a potential source of systematic error. For an 
asymmetric beam, one could make the optimistic assump¬ 
tion beam multipoles in a given frequency channel could be 
written as the product of a frequency-dependent, symmetric 
part fepo)(r) which is assumed to be the same for all detec¬ 
tors in the band, and a frequency-independent, asymmetric 
part which may depend on the detector: 

bd(im)(y) ~ Z>(zo)(T)[^d(z m) M/Z<qo)M], (31) 
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where uo is the central frequency in the band. We could 
then solve for the multipoles of the sky convolved with an 
effective symmetric beam with multipoles bfi m ){n) and inte¬ 
grated across the frequency band with filter v(n). By factor¬ 
ing the beam in this way, we can simultaneously improve the 
condition number of A^N X A and remove most of the bias 
due to variation of the main beam width with frequency. 


3 SIMPLE SCAN STRATEGIES 


In this section we investigate the structure of the inverse 
covariance matrix, C -1 , of the recovered multipoles for two 
simple scanning strategies. The first is a random pointing of 
the spin axis for each ring, which ensures uniform coverage 
of the entire sky. The second is constant latitude scanning, 
where all rings have the same polar angle, 9 r = 9. Constant 
latitude scanning is a useful approximation to the scan strat¬ 
egy of the Planck satellite, where the spin axis stays close 
to the plane of the ecliptic. By adopting coarse restrictions 
on the behaviour of the instrument, we are able to demon¬ 
strate that, under these restrictions, the harmonic estimate 
of the sky, equation (28), is statistically equivalent to that 
obtained via conventional map-making and spherical anal¬ 
ysis. By relaxing some restrictions on the behaviour of the 
instrument, we obtain some useful extensions of the known 
results. 

Throughout this section we impose a number of limits 
on the instrument response: (i) The spin axis position is al¬ 
ways aligned with the z-axis of the instrument reference sys¬ 
tem (see Section 2.2), so that the ring opening angles, a r d, 
are independent of r, and the focal plane rotations K r d = 0 ; 
(ii) There is no variation in average spin rate between rings, 
so we can drop the subscript from u> r ; (iii) The band width 
of the spectral filters is sufficiently narrow that 

J a(iL)nb(tm')(^)rH)du = af l ^ l) (ty 0 )b[ lrn , ) {i' o) 

x f v(n)dn, (32) 


where no is the central frequency of the band; (iv) Bolometer 
time constants, Td, and the sampling period, A d , are such 
that TdtOrhnax <S 1 and Adutrlmax <S 1 , in which case the 
effect of the instrument response can be approximated by 
simply a gain at zero lag; (v) The noise power spectrum is 
white so that 


N( r dn)(r'd'n') — S rr ' 5dd' &d ; 


(33) 


where <7% is a constant for each detector, which can be related 
to the detector sensitivity, Sd With the restrictions above, 

r i—i 


Sd 


A e ff,dGdT 


dB(v 0 ,T) 


dT 


v(y ) &v 


^ The detector sensitivity is defined such that for an integration 
time At with a single pointing, the excess signal to noise is unity 
when the instrument is illuminated with an unpolarized black- 
body sky with a temperature excess Ts/VAt over the average 
CMB temperature, T. If the units of the gain Gd are e.g. VW _1 , 
so that Srdi'ift) is a readout Voltage, ad will also have the dimen¬ 
sions of a Voltage, and Sd will have units of e.g. (p.K/K)s 1//2 . 


X 



(34) 


where B(n,T) is the Planck function and the average CMB 
temperature is T. Note that a consequence of restriction (i) 
is that Bf d(ln) is independent of ring, Bf d(ln) = B d(Jn) . If we 
define dimensionless multipoles by 

-l 

o), (35) 

the inverse of the dimensionless error covariance matrix, 
C , has components 


l(Zm) 


f dB(v 0 ,T) 


dT 


C 


-1 

(PZm)(P'Z'm / ) 


^ 5 ] S ~d 2 { Dl mn{4>r, 9 r , 0)B% ln) 

rdn 

X [D rn > n {<j>r,0r,0)Bd(l’n)] }i (36) 


where the sum is over detectors, rings, and |n| ^ min (1,1')- 
We now proceed to analyse C 1 for randomly-positioned 
rings and constant latitude scanning. 


3.1 Randomly-positioned rings 


We consider the limit of equation (36) as the number of 
rings N r —> oo while keeping the survey length finite, in 
which case we can replace the sum over rings by an appro¬ 
priate integral. In practice, for fi nite N r th e continuum limit 
should hold for scales with l ^ HN T /( 47 r). Placing the rings 
at random is equivalent to enforcing uniform coverage of the 
full sky, which will allow us to make contact with existing an¬ 
alytic results obtained from the map-based formalism. The 
sum over randomly-positioned rings in equation (36) can 
be replaced by [A r /( 47 r)] J dfl r , where dfl r = d<(> r d cos # r . 
Making use of the orthogonality of the D l mm , (a, (3, 7 ) over 
the SO(3) group manifold (Brink & Satchler 1993), we find 


A i ,in(')>r-,6 , r,0)[Rj n ,„(<(> T .,6» r ,0)]* dn r = 5 W 5 mm ’ ,(37) 


so that C 1 simplifies to 

C(plm,)(P'l'm’) = bw^mm 1 .y ~ J s d &d(ln)B d(;„) , (38) 


where the total mission time is T m = 2-kN s N t / ui. Equa¬ 
tion (38) can be further simplified by using equation (13) 
to substitute for B d (i n \ and then performing the sum over n 
noting that [D^^a, /?, 7 )]* = D l nrn {- 7 , -/?, -a). The result 
is 

/'■ 1 _ xx T m 

~ 21 + 1 

^ ^ ' N/ b d (im") (yofbd(im") (^o)- (39) 

dm." 


Note that there is now no dependence on the opening angles 
of the detectors, a r d = ay • This is a consequence of assuming 
white noise. If the noise had a characteristic timescale, this 
would combine with the spin velocity and opening angle to 
define a characteristic angular scale for the projected noise, 
and C would then depend on ad- 

To make further progress, we consider initially the case 
where the beams are axisymmetric, and have no cross-polar 
contamination. In this limit, the intensity multipoles can be 
written in terms of a window function, Wi,d- 
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t>d(lm){ v o) ~ 



Wl’dSmO, 


and the linear polarization multipoles in terms of spin- 
weight 2 window functions, 2 Wi, d , (Challinor et al. 2000; 
see also Ng & Liu 1999): 


bd(lm){ u 0 ) — 


)M = 


l^-2W l , d e~ impd ( 

oZtt 


5m2 Sm. — 2 ), (41) 


-iJ^^2 w i, d e- impd (5 m 2-Sm-2). (42) 


Here, p d is the angle between the polarization direction on 
axis and the normal to the plane containing the main beam 
direction and the spin axis. (The unit normal to this plane is 
<r y when the detector is in its horn’s reference configuration). 
By assumption, the circular polarization multipoles for the 
beam vanish. Substituting into equation (39), we find 




= Sw5, 


ll'Omm'Opp ' 


X W d \Sf Wl d + (Sq + Sc)2Wl t d], (43) 


where we have introduced the weight per solid angle (Knox 
1995), which for uniform coverage of the sky is w d = 
Tm/(47rs|). Note that the covariance matrix is diagonal, and 
that it does not depend on the relative orientations of the 
polarization directions of the various detectors. This is not 
surprising since for randomly-positioned rings, every point 
on the sky is traversed by every detector in every possible 
orientation, and we have assumed the instrument noise to 
be uncorrelated between detectors. 

It is straightforward to show that the covariance matrix 
in equation (43) is equal to that obtained by constructing 
optimal (minimum variance), beam-smoothed maps of the 
Stokes parameters for each detector, and then optimally ex¬ 
tracting the from a joint analysis of the maps. If the 
main beam of the dth detector lies in the pth pixel when 
the ring phase is ip on the rth ring, the signal contribution 
to the phase-ordered data, equation (17), can be written as 
(see e.g. Challinor et al. 2000) 


^rd(V’) — -4- e ff jd G d T 


dB(u 0 ,T) 


v(v) dp[/eff,d(e p ) 


- Qett,d(e p ) cos 2?? + U e ff, d (e p ) sin 2r/]. (44) 

under the restrictions described above. In equation (44), e p 
is in the direction of the pth pixel, and has polar angle 9 P 
and azimuth (f> p , and the angle p is the angle between the 
polarization direction (on the sky) along the beam axis and 
the plane containing e p and the 2 -axis of the fixed reference 
system. The beam-smoothed fields are 

7eft,d(e) = ^2 Wi, d ali m )Y( lm )(e) (45) 

lm 

-^={Q e a,d ± iU e s,d)( e ) = ^[ 2 Wi,d(afi m ) Tiafim)) 

lm 

X y 2 V( !m )(e)], (46) 

where the ±2 Y(im) are examples of the spin-weight harmon¬ 
ics, defined for integer s by (Goldberg et al. 1967) 


D l m _ s (a,(3, 7 ) = {-l) s J^T[sY( lrn) {P,ay- 


For randomly-positioned rings in the limit N r —» 00 , the 
angles ip are uniformly covered in any given pixel. Given 
our assumptions, the optimal, unbiased estimate for I e g,d in 
the pth pixel is given by direct averaging of the appropri¬ 
ate t r d(ip)- For Qeff.d the t rd [ip) are averaged with weight 
—2 cos 2p, and for U e g, d the weight is 2 sin 277 . The errors 
on the smoothed maps are uncorrelated between Stokes pa¬ 
rameters, and between pixels. The errors on these maps have 
weights per solid angle w d for I e s.d, and w d /2 for Q e ff,d and 

U e H,d- 

For uncorrelated noise between Stokes parameters, pix¬ 
els, and detectors, the maximum-likelihood estimate of the 
®(im) f rom H le smoothed maps reduces to minimising the 
(absolute) squared residuals between the left and right-hand 
sides of equations (45) and (46), with each pixel for each de¬ 
tector entering with the appropriate weight per solid angle. 
In the case of uniform coverage considered in this subsection, 
the maximum-likelihood estimates of the d^ lm ^ are equiva¬ 
lent to performing a spherical transform of the smoothed 
maps and then averaging these across the detectors, with 
each detector carrying statistical weight Wi )d w d . However, 
it will be useful for the next subsection to note the form that 
the error covariance matrix takes when we relax the assump¬ 
tion of uniform coverage, while retaining the assumptions 
of uncorrelated errors between Stokes parameters, detectors 
and pixels. The weights per solid angle for I e s,d are then 
pixel-dependent, w d , P , and the weights for Qeff.d and U e g, d 
are u>d,p/2. The non-vanishing entries of the inverse of the 
(Hermitian) error covariance matrix are then: 

C(llm)(IVm') = ^ ^ Wl,dWip d Wd,A (lm)Y(Vm') Aflp, (48) 

pd, 

C(Glm)(GVm') = 'y\2W h d2W V ,dWd, P \{2Yf lm) 2\\ Vrn , ) 

pd 

+ -2Yf l m)-2Y(l, rn >))AQ p \, (49) 

C(Glm)(CVm') = iy\2Wl,d2W V , d Wd, P ^YC lm) 2Y^ ml) 

pd 

— - 2 K(; m )_ 2 Y(i' m '))AfIp]j (50) 

and C(Clm)(Gl'm') = C(Zlm)(Gl'm'Y Here > is the S ° lid 

angle subtended by the pth pixel. For the case of uniform 
coverage, w d , P = w d , the orthonormality of the spin-weight 
harmonics reduces equations (48)-(50) to the result derived 
from the harmonic model, equation (43), if we ignore pix- 
elisation effects. For this simple example, the estimates of 
the from the harmonic model and the map-making 

route are both unbiased and have the same error covariance. 
It follows that they are statistically equivalent at the level 
of second-order correlations (and at all orders for Gaussian 
noise). 

3.1.1 Effect of beam asymmetry 

We now consider the effect of a known beam asymmetry 
on the covariance matrices obtained via the harmonic and 
map-based routes. For simplicity, we shall only consider un¬ 
polarized detectors with bivariate Gaussian profiles, with ec¬ 
centricity e d - Consider initially a detector in its horn’s refer- 
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ence configuration, with the major axis of the iso-directivity 
ellipse aligned with cr x . For beam widths ad -C 1, we have 

(51) 

which has non-zero multipoles, in the limit of large l, 


Id{e\ vo) 


exp{—0 2 [cos 2 0 + sin 2 0/(1 


2naj^/T^e 


bd(lm) 


2^ + 1 r / 7 2 2 2 j - 

tm/ 2 \J' a e d /4)e 


l 3 ^(l-e 2 d /2)/2 


(52) 


for m even, where 7 m (x) are modified Bessel functions. Note 
that beam asymmetry is only significant for l 1 /{ade d )- 
The multipoles for a detector whose iso-directivity ellipse 
has major axis at angle -y d to a x pick up an additional phase 
exp(— im-yd)- 

The expression for (equation 39) involves 

the sum ^/\, <( l^d(in)( t 'o)| 2 - For large l , this is easily evalu¬ 
ated by substituting the integral representation of the modi¬ 
fied Bessel functions (or, equivalently, azimuthally averaging 
the absolute square of the flat-space Fourier transform of the 
directivity). The result is 


12 2Z + 1 r n2 _ 2 _2 /nN -l 2 <7+1-e 2 ./2) 


Y, \bd(i n ){v 0 )\ 2 = --—i 0 (? 2 o-Je 2 /2)e 


(53) 


\n\^l 

so that 
C~ 1 

*■"' (Ilm)(Il'm') 


Sll'Smm 1 "y w d Io{l 2 ade 2 d /2) 


x e 


(54) 


If a known beam asymmetry is to be accounted for with 
pixel-based map-making techniques, it is necessary to in¬ 
clude the asymmetry in the real-space pointing matrix. A 
simpler procedure suggests itself for the case of randomly 
positioned rings, where we can exploit the fact that every 
pixel is sampled with each detector in every orientation. Av¬ 
eraging the data in each pixel from a given detector gives an 
unbiased estimate of the pixelised sky smoothed with an ef¬ 
fective window function. This window function follows from 
azimuthally averaging the beam, so that 


Wea,i,d = yj = /o(Z 2 ^el/4)e- i2 ^ (1 -^ /2)/2 .(55) 

The analysis of the smoothed maps proceeds as in Sec¬ 
tion 3.1, so that the inverse covariance matrix evaluates to 

C(nm)(ii'm') = bivSmm'y^Wdllifalel/A) 

d 


x e 


-J 2 <r 2 (l-e 2 /2) 


(56) 


Clearly there is some information loss on averaging the data 
in the manner described. For a single detector, the ratios of 
the diagonal elements of the covariance matrices obtained 
from equations (54) and (56) are 


M^e 2 /2) 


= 1- 7^d°de d ) A + ■ ■ ■ ■ 


(57) 


biased estimate of the sky in the presence of a beam asym¬ 
metry. This will be most acute for scanning strategies where 
a large number of pixels are sampled with the detectors in 
only a narrow range of orientations. Such strategies include 
constant latitude scanning, to which we now turn. 


3.2 Constant latitude rings 

To analyse constant latitude scanning, we consider equa¬ 
tion (36) under the restriction 6 r = 9. Following the discus¬ 
sion at the start of Section 3.1, we consider the continuum 
limit. For constant latitude scans we replace the sum over 
rings by [7V r /(27r)] f d0 r . Expanding the D-matrices as in 
equation (9) and performing the integral over <j> r , the ex¬ 
pression for C 1 reduces to 

67 Plm)(P'l'm') 47T(5 mm / ^ ^ W d d rnn (0)d rn / Jl jO)Bd{ln) Bd(l'n) ? 

dn 

(58) 

where the sum is over detectors and \n\ ^ min (1,1'). We 
obtain uncorrelated errors between modes with different m 
due to the azimuthal symmetry of the sky coverage (Oh et 
al. 1999). Since C _1 is block-diagonal, it can be inverted in 
O(Zmax) operations. I 11 practice, variations in the instrument 
parameters through the mission will spoil this exact block- 
diagonal structure, as will any precession in the latitude 9. 
However, approximating C _1 by its diagonal blocks may still 
provide an adequate preconditioner for e.g. a conjugate gra¬ 
dient reconstruction of the sky. 

In equation (58), w d = T m /(4-7rs 2 ) refers to the weight 
per solid angle for uniform coverage of the whole sky, which 
differs from the pixel-dependent w d , P due to the variation in 
time spent per solid angle. By differentiating the geometric 
relation 


cos 9 P = cos ad cos 9 — sin ad sin 9 cos ip, 


(59) 


which relates the polar angle 9 P of the pixel containing the 
main beam to the ring phase ip, it is straightforward to show 
that 


2 

Uld,p — 

7T 


®{B d , P ) 

\JBdf> 


w d , 


(60) 


where 0(x) is the Heaviside unit step function, and we have 
introduced the quantity 

Bd, P = 1 —cos 2 ay — cos 2 9 P — cos 2 9 + 2 cos ad cos 9 V cos 0(61) 


for later convenience. 

To proceed further, we make the assumption that the 
detector sensitivities and ring opening angles are all equal, 
so that w d = w, ad = a, and B d:P = B p . It follows that the 
pixel-dependent weight per solid angle is also independent of 
detector: w d , P = w p . The sum over detectors in equation (58) 
then reduces to 

^ ^ BdO + Bda' n) f (^)^nm'" (®) 

d m"m" r 

X bd(lm")bd(l*m"')\- (62) 

d 


This information loss is only significant on scales below that 
of the beam asymmetry (« a d e d )- For more general scan¬ 
ning, averaging data in a pixel will necessarily produce a 


We further assume that all detectors have axisymmetric, 
co-polar beams which are equivalent up to rotations about 
the beam axis. In this case, the multipoles are given by 
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equations (40)-(42) with the angles pd determining the 
relative polarization orientations. The choice of angles pd 
has a strong impact on the covariance structure of the 
estimated multipoles. Designing the focal plane so that 
exp(2 ipd) = 0 ensures that the errors are uncorrelated 
between the total intensity and the linear polarization. The 
correlation between the G and C multipoles is controlled 
in part by the sum y d exp(4 ipd), although demanding that 
this sum vanish is not a sufficient condition to ensure un¬ 
correlated errors between G and C. However, the condition 
exp(4 ipd) = 0 does ensure uncorrelated errors between 

the smoothed Stokes fieldsll Q e ff(e) and U e fr(e) which, un¬ 
der the stringent restrictions adopted in this section, can be 
estimated in an optimal manner by a \ 2 fitting of the data 
from all detectors obtained when the main beam of each 
lies in the given pixel (Couchot et al. 1999). The smoothed 
maps estimated in this manner have uncorrelated errors be¬ 
tween pixels with weights per solid angle NdW p for 7 e fj, and 
NdWp/2 for Qefi and U e s, where, recall, Nd is the num¬ 
ber of detectors. Following Couchot et al. (1999), we de¬ 
note arrangements of polarimeters with y exp(2 ipd) = 0 
and y^ exp(4 ipd) = 0 as optimal configurations. In such a 
configuration, the non-vanishing components of the inverse 
covariance matrix evaluate to 

= Smm’WNdWlWl’ \J(21 + 1)(2Z' + 1) 

X y: d l mn (e)d l mln {0)d l nO (a)d l nO (a ), (63) 

n 

C(imG)(Vm'G) = Smm'WNd 2 Wi 2 Wy y/(21 + 1)(2F + 1) 

X 5 ~2{d l mn{0)d l m ' n {9)\ [d l n2 ( a)d l n2 (a) 

n 

+ d! n - 2 (a)d l n - 2 (a)]}, (64) 

C(lmG)(l'm'C) = —ifimm'WNd 2 W 1 2 Wp \J(21 + 1 )( 2 V + 1 ) 

X y^Xd l mn {e)d l mln (e)\[d l n2 (a)d l n2 {a) 

n 

- d! n _ 2 {a)d l n _ 2 (a)]}, (65) 

with CfTmCXi'm'C) = ^(TmGXi'm'G)- In equations (63)-(65), 
the sums are over |n| ^ min (1,1'). 

We show in Appendix B that the sum over n on the 
right-hand sides of equations (63)-(65) can be reduced to 
matrix elements of the time spent per solid angle with the 
appropriate spin weight basis: 

fimm 1 ^ ' \J (2 1 + 1)(2 1' + 1 )d l mrl (8)d l rn r n (9)d l ns (a)d l na (a) 


2 Q(B P 

7T 




m ) (^p) — A(I'm') (^p) dflp, 


( 66 ) 


for integer s. Using this result in equations (63)-(65), and 
recalling equation (60), the components of the inverse co- 


II The smoothed fields are defined in equations (45) and (46). 
We have dropped the subscript d since we are assuming that the 
beam window functions are the same for all detectors. 


variance matrix can be reduced to 
G~ 1 

^ (Ilm)(Il r m') 

n~ i 

^ (Glm)(GV m') 


n- 1 

Glm)(Cl'm') 

— -2U(im)-2T(i'm'))dfip]. (69) 

Note that these results automatically take account of in¬ 
complete sky coverage, through the presence of Q(B P ). The 
function B p < 0 in pixels that are not sampled by the main 
beam during the mission, so such pixels make no contribu¬ 
tion to the covariance matrix. As noted in Section 2.5, the 
presence of unsampled regions can render CD 1 numerically 
singular. For the Planck mission the nominal ring opening 
angle a « 85°, which would leave 5° holes at the ecliptic 
poles for constant latitude scanning in the ecliptic. How¬ 
ever, a proposal to precess the spin axis out of the ecliptic 
plane with an amplitude of 10° has the benefit of providing 
complete sky coverage, which would render CU 1 invertible 
(although no longer block-diagonal). 

To compare equations (67)-(69) with the results ob¬ 
tained from a pixelised map, we can make use of equa¬ 
tions (48)-(50) which, it will be recalled, give the errors on 
the maximum-likelihood solution for the multipoles obtained 
from Nd sets of smoothed maps. Here, we have only one 
map (estimated using data from all detectors as described 
above), with weights per pixel NdW p for J e g and NdW p /2 
for Qeff and U e a- With these conditions, it is straightfor¬ 
ward to verify that equations (48)-(50) reduce to pixelised 
versions of equations (67)-(69). This observation confirms 
the statistical equivalence of the map-based and harmonic 
routes through to the multipoles of the sky for constant lat¬ 
itude scanning, under the strong restrictions adopted in this 
section. 

More generally, since the maps and their multipoles con¬ 
tain the same information (disregarding pixelisation errors), 
it will always be the case that the optimal map-making and 
map-to-multipole algorithms are statistically equivalent to 
the optimal harmonic estimate of the multipoles. However, 
as demonstrated in Section 2, the harmonic model provides 
a more natural framework for the inclusion of a number 
of systematic effects in the instrument modelling and data 
analysis pipeline. The inclusion of some such effects is es¬ 
sential to maintain the integrity of the final data products. 


WiW v J w p Y'(im)Y(i'm')<ifl p , (67) 
2 Wl 2 Wl' j [w p i(2F(L)2lW) 

+ -A(lm)- 2 Y(Vm')) dfip], (68) 

i 2 Wi 2 Wi> f[w p ^( 2 Y^ m ^ 2 Y(ifm') 


4 DISCUSSION 

In this section we discuss a number of important issues that 
arise during the map-making stage. The treatment of low 
frequency noise (‘destriping’) in the harmonic model is de¬ 
scribed in Section 4.1. In Section 4.2 we comment on the 
control of scan-synchronous instrument effects, focusing on 
strayliglit from bright sources picked up through the side- 
lobes of the telescope, and the modulation of the dipole in 
the rest-frame of the experiment. Finally, the overall cali¬ 
bration of the experiment to external standards is described 
in Section 4.3. 
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4.1 Destriping 

In Section 2.4 we established that for large N s , the main im¬ 
pact of noise power at frequencies below the spin frequency 
is concentrated at the n = 0 Fourier modes of the phase- 
ordered data (see also, e.g. Delabrouillc 1998). If the noise 
power varies sufficiently rapidly in the range 0 ^ ui ^ u r /N s 
these offsets may be correlated between rings. Careful treat¬ 
ment of low frequency noise is essential to avoid undesirable 
long range noise correlations (‘stripes’) in the maps. Given 
the difficulty of inverting the time-time noise covariance 
matrix in the presence of low frequency noise, additional 
pre-processing steps are usually performed prior to map¬ 
making. Pre-whitening the time- or phase-ordered data (e.g. 
Tegmark 1997) uses a prior knowledge of the noise power 
spectrum, A/d(w), to represent the data in a basis where 
there are no noise correlations. In essence, this involves the 
subtraction of an optimal offset from each ring of data. Na- 
toli et al. (2001) have successfully applied this method to the 
30 GHz channel of the Planck Low Frequency Instrument 
under the assumption of a symmetric beam profile. Other 
destriping methods advocated for the Planck mission, e.g. 
Delabrouille (1998) and Revenu et al. (2000), also involve 
the subtraction of offsets from the rings, but these offsets 
are obtained directly from the data by exploiting the inter¬ 
sections of rings. 

The harmonic model represents the phase-ordered data 
in the Fourier basis on rings which ensures that the noise 
covariance matrix, l\l, is very sparse, even in the presence of 
significant low frequency power. This removes the need for 
any additional pre-whitening. A potential problem arises if 
the noise power spectrum is not known accurately at low 
frequencies, since this may compromise the quality of the 
sky reconstruction. The harmonic model suggests a simple 
solution to this problem: remove the n = 0 Fourier modes 
from the analysis. From a Bayesian viewpoint, we must now 
infer the df lm ^ from a knowledge of t( r dn) with n > 0. This 
inversion requires the probability density function (pdf) for 
n( r dn) with n > 0, which is obtained from the full pdf by 
integrating over the n = 0 modes. However, for Gaussian 
noise, the pdf factors into products of pdfs for the individ¬ 
ual Fourier modes as a consequence of equation (24), so that 
integrating over the n = 0 modes is trivial. Note also that 
since the n = 0 modes of the data are the only ones that 
depend on the l = 0 modes of the total intensity and circular 
polarization, these monopole modes can no longer be recon¬ 
structed. To implement the scheme one can either reformu¬ 
late the problem with the n = 0 data modes removed at the 
outset, or, equivalently, let Af( r rf 0 )(r- / iio) —> oo, and perform 
the inversion in the subspace orthogonal to a[ 00 ) and d(g 0 y 
Removing the n = 0 modes does not bias the solution, but it 
is no longer optimal given all the data, since t( r do) receives a 
contribution from all multipoles. However, for large l mSLX the 
information loss should be small for the l > 0 multipoles. Ex¬ 
cluding the n = 0 modes provides a very robust, simple way 
of removing undesirable long range correlations in the maps 
that might arise if the low frequency noise were estimated 
inaccurately. Note also that this ‘destriping’ method handles 
the polarization signal seamlessly, avoiding the technicalities 
involved in estimating offsets from ring intersections for po¬ 
larized detectors (Revenu et al. 2000). Removing the n = 0 
and l = 0 modes from the analysis also neatly solves the 


problem of degeneracy that arises if the experiment is in¬ 
sensitive to the monopole - this situation may arise for the 
Planck HFI since the favoured readout electronic system is 
insensitive to very low frequencies (Gaertner et al. 1997). 

4.2 Control of Scan-Synchronous Instrument 
Effects 

Instrument effects that produce signals synchronous with 
the rotation of the instrument will not be suppressed by 
co-adding the time-ordered data to form the phase-ordered 
data. Given complete knowledge of the instrument response, 
these potential sources of systematic error can be controlled 
by adopting a more refined model of the instrument response 
during the reconstruction of the sky. In this subsection we 
discuss the removal of two potential systematic effects within 
the context of the harmonic data model. 

4-2.1 Sidelobe Reconstruction 

An important source of systematic error is straylight from 
e.g. the Galaxy and the CMB dipole entering the instru¬ 
ment through the sidelobes of the telescope. Here, the dif¬ 
ficulty lies not in the inclusion of the sidelobes in the re¬ 
construction process, since the harmonic model allows for a 
complete description of the beams, but rather in the lim¬ 
ited knowledge of the sidelobes that will be available from 
simulations and from in-flight calibrations (Burigana et al. 
2001; van Leeuwen et al. 2001). Although the sidelobes can 
be expected to be highly polarized, the direction should be 
sufficiently random that it is adequate to model only the 
directivity Id(e;u o) in the sidelobes. 

Delabrouillc et al. (1998b) have proposed an iterative 
scheme for estimating sidelobe corrections. In their method, 
an estimate of the sky is obtained from a first pass of the 
map-making algorithm ignoring sidelobe corrections, and 
this is then used to compute the difference between the 
observed time-ordered data and the contribution of the 
sky coming through the main beam. This difference is at¬ 
tributable to instrument noise and the sidelobe signal, and 
can be inverted (with suitable regularisation) to estimate 
the values of the directivity in the sidelobe pixels. With this 
improved knowledge of the beam, and the original estimate 
of the sky, an estimate of the sidelobe contribution can be 
subtracted from the time-ordered data, and an improved es¬ 
timate of the sky obtained from a further pass of the original 
map-making algorithm (again with sidelobes ignored). This 
process can be iterated to give consistent estimates of the 
sky and the directivity in the sidelobe. It should be possi¬ 
ble to implement a similar scheme in the harmonic model, 
with the sidelobes parameterised by a modest number of 
multipoles and the sidelobe contribution estimated from the 
t( r dn) i although we have not attempted this yet. One poten¬ 
tial problem with such a scheme is that the sidelobes may 
contain rather sharp features due to geometrical effects, in 
which case a parameterisation of the sidelobes with spherical 
multipoles may not be ideal. The monopole 6( 00 j (no) should 
not be varied during sidelobe reconstruction since its value 
is fixed at 1 /by definition. The iterative reconstruction 
of the sidelobes does not require an absolute calibration of 
the detector gains, and should be performed prior to the 
calibration procedures described in Section 4.3. 
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4-2.2 Dipole variation 


One subtlety that we have overlooked so far, which is sig¬ 
nificant for survey missions from space, is the variation of 
satellites (linear) velocity during the year. The orbital ve¬ 
locity of the earth about the sun ~ 30kms -1 , which induces 
a yearly modulation in the multipoles seen in the satellite’s 
rest frame. For sensitivities equivalent to the Planck mis¬ 
sion (a few fiK/K in lOOarcmin 2 pixels at 100 GHz, for 
twelve months of observation), the most significant effect is 
the modulation of the intensity dipole due to the monopole 
with amplitude ~ 10 -4 K. If we denote the multipoles in the 
frame of the satellite by B (z/), and continue to denote 
the multipoles on the { a x , <r y , a z j basis fixed relative to the 
solar system by ap m )(^), we have 


afim), eM ~ (pim)^) + df <5fi 


yind v 

x l V 


a (oo) ( u ) 


(70) 


where v is the (time-dependent) velocity of the earth rela¬ 
tive to the sun, and e is a unit vector. The B (z/) are 
defined relative to a basis obtained by Lorentz boosting the 
{ <T x ,a y ,cr ,} basis. Aberration effects on the high l multi¬ 
poles may also be significant for Planck (Challinor & van 
Leeuwen, in preparation). 

It is the time-dependent multipoles ag, m ) E (obtained by 
averaging B (z/) across the spectral filter) which should 
now appear in equation (20), relating the signal to the sky. 
However, we must still solve for the time-independent mul¬ 
tipoles • If the monopole a 4 00 ) (u) were available at the 
start of the analysis for all frequency channels, the contribu¬ 
tion of the time-dependent part of aq m ) e(z') to the signal 
could be subtracted from the t( r dn) prior to solving for the 
multipoles. Note that this process only affects the n = 0 
and n — 1 Fourier modes, t( r d n ), but these influence all 
multipoles of the reconstructed sky. More realistically, the 
monopole a^ 00 ^(iy) may only be poorly determined prior to 
the analysis, in which case it will be necessary to iterate the 
estimation of the multipoles, using the monopole determined 
in the previous iteration to correct the t( r dn). This process 
is clearly not possible if the n = 0 modes are ignored in the 
analysis, as may be necessary for the reasons discussed in 
Section 4.1. A radical way to circumvent the problems of 
poorly determined low frequency noise and dipole modula¬ 
tion is to remove the n = 1 modes from the analysis also. 
However, an internally reconstructed dipole appears to be 
essential for establishing an absolute normalisation of the 
experiment (see next subsection). 


4.3 Gain calibration 

The methods described by van Leeuwen et al. (2001) for 
producing internally consistent ring-sets do not provide an 
absolute gain calibration for each detector. By forgoing the 
introduction of any external flux calibrator during the re¬ 
construction of the ring-sets, only the ratios of the gains of 
all detectors in a given frequency band will be known; the 
absolute gains will only be determined up to an overall fac¬ 
tor, g. This factor is defined to be the ratio of the true gains, 
Gd, true, to those assumed in the model of the instrument, Gd 


(e.g. equation 19). An external calibrator provides prior in¬ 
formation on the sky, Pr(a), which can be used to constrain 
g. For simplicity, we assume a Gaussian prior: 

Pr(a) oc |7rC ex tr 1/2 exp[-(a - a ex t) t CJ x 1 t (a - a ex t)/2],(71) 


where a ex t is the most likely prior sky, and C ex t is the as¬ 
sociated error. The sky and normalisation can be now es¬ 
timated directly from the (Fourier) ring data, t , by max¬ 
imising Pr(tja,g)Pr(a) (assuming a uniform prior on g). In 
the usual limit where the experiment determines the sky to 
much better precision than the calibrator, the best estimate 
of the sky is a/g, where a is the maximum-likelihood esti¬ 
mate for the sky given in equation (28), obtained with gains 
Gd, and the best estimate of the normalisation is given by 


9 = 




®(®^C ext a ex t) 


(72) 


In practice, C/^aext and C/j/t may have to be estimated from 
existing observations on patches of the sky. In such cases, 
equation (72) reduces to the solution obtained by minimis¬ 
ing the (error-weighted) residuals between the reconstructed 
map, X/Zm 9~ 1 ®(im)Y(im){ e ), and the calibration map. The 
requirement that the calibration maps be at the same fre¬ 
quencies as the instrument channels can be removed by se¬ 
lecting a large patch, at high Galactic latitude, where the 
large-scale signal is dominated by the CMB dipole. 


5 FURTHER ANALYSIS OF FREQUENCY 
MAPS 

Although the main focus of this paper is the reconstruc¬ 
tion of frequency maps of the sky in multipole space, in this 
section we offer some comments on the subsequent analy¬ 
sis of these maps into their astrophysical components, and 
estimated power spectra. 


5.1 Component separation 

The frequency maps contain contributions from a number 
of astrophysical components which must be separated on 
the basis of their assumed frequency spectra (and possibly 
power spectra). Unresolved radio sources require separate 
processing since they cannot be accurately modelled as a 
population with a single frequency spectrum. If bright point 
sources were not removed from the ring-sets prior to map¬ 
making, their contribution should be filtered from the mul¬ 
tipoles before attempting component separation (Vielva et 
al. 2001). Within the timescale of Planck, the majority of 
point sources expected to be visible in the 857 GHz chan¬ 
nel of the HFI should already have been catalogued prior to 
launch by the ASTRO-F survey **. Such catalogues will be 
valuable for the geometrical calibration of the ring-sets using 
the positions of identified point-sources (van Leeuwen et al. 
2001), and also for the removal of point sources from the re¬ 
constructed frequency maps. If bright point sources have al¬ 
ready been removed from the ring-sets prior to map-making, 
component separation can proceed directly from the recon¬ 
structed multipoles. However, statistical reconstruction of 


** http://www.ir.isas.ac.jp/ASTRO-F/ 
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faint point sources may still be desirable, in which case the 
faint point sources can be included as a generalised noise in 
the separation algorithm for the other astrophysical compo¬ 
nents (Hobson et al. 1999). 

The full-sky, maximum-entropy component separation 
algorithm developed recently by Stolyarov et al. (2001) per¬ 
forms the separation in the spherical multipole basis, so that 
it can use the outputs of the multipole estimation process 
described here directly. Although the algorithm of Hobson et 
al. does not handle polarization in its current form, the ex¬ 
tension to polarized components should be straightforward. 
(The Wiener separation of polarized components has been 
implemented successfully by Bouchet, Prunet & Sethi 1999 
for small patches of sky.) A more demanding problem for 
high resolution, all-sky separation algorithms is the inclu¬ 
sion of non-isotropic noise, since then the separation cannot 
be performed multipole by multipole due to the correlated 
errors (e.g. Prunet et al. 2001). 


5.2 Power spectrum estimation 


Power spectrum estimation can also be performed efficiently 
in multipole space. The simplest, unbiased estimator which 
is quadratic in the estimated dimensionless multipoles &(; m ) 
is of the form 


cr' = 


C 21 + 1 ) 




(73) 


The variance of this estimator is calculated in Kamionkowski 
et al. (1997) under the assumptions of full sky coverage and 
isotropic pixel noise which is uncorrelated between Stokes 
parameters. In this case, the quadratic estimator in equa¬ 
tion (73) is equivalent to the maximum-likelihood estimator. 

For more general noise properties and scan strategies 
the estimator in equation (73) is sub-optimal. Maximum- 
likelihood techniques have been developed which allow ac¬ 
curate computation of the temperature power spectrum in 
0{ltia,x) operations for scanning strategies with (near) az¬ 
imuthal symmetry (Oh et al. 1999). This method cannot be 
directly in the context of the harmonic model since Oh et al. 
store the inverse (noise) covariance matrix, C -1 , in the pixel 
basis only, where they assume it is diagonal. Forward and 
inverse fast spherical transforms are then employed to apply 
C 1 efficiently in multipole space. Applying C _1 directly in 
multipole space would not increase the overall operations 
count significantly, but would increase the memory require¬ 
ments to O(^max) hi general. 

One potential issue in power spectrum estimation is the 
treatment of regions near to the Galactic plane. Although 
component separation algorithms can reconstruct the CMB 
very well even at low Galactic latitude (Stolyarov et al. 
2001), it may still be necessary to remove traces of the 
Galaxy by brute force. Given a pixelised map, the pixels near 
the plane can be excised from the analysis before estimating 
the multipoles and their error covariance, and subsequent 
power spectrum estimation. (Alternatively, power spectrum 
estimation can be performed directly in real space away from 
the Galactic plane, although the lack of a sparse signal co- 
variance matrix in the pixel representation is problematic 
for high resolution data-sets.) 

Within the context of the harmonic model, the Galac¬ 
tic plane can be dealt with by projecting the reconstructed 


multipoles, dq m ), into a subspace which is (almost) free 
from Galactic contamination. For simplicity, we shall only 
describe the procedure for a total power measurement; the 
extension to polarized data is straightforward. We define the 
Hermitian matrix P to have components 

= I Y(lm) (^)^ (Vm') (^) dfl, (74) 

Js 2 ' 

where the integral is over the region of the sky that we 
wish to retain in the analysis. In the limit I max —> oo, 
P(im.)(i'm') are the matrix elements of a projection operator 
(i.e. P 2 = P) which projects functions into the region S 2 '. 
For finite Z max 1. P is almost a projection operator since 
the distribution of its eigenvalues is almost bimodal with val¬ 
ues clustered close to zero and one (Mortlock et al. 2001). 
The fraction of eigenvalues close to unity is approximately 
the fraction of the sky retained; the remainder are nearly 
all zero. Performing (maximum-likelihood) power spectrum 
estimation with the data object Pa would not remove the 
Galactic contamination since P is (formally) invertible. To 
enforce rejection of the Galaxy, we introduce a projection 
operator P which is obtained from P by retaining the eigen¬ 
vectors but setting those eigenvalues greater then 1 — e to 
unity, and to zero otherwise. The choice of threshold, e ^ 1, 
must be determined from simulations to ensure good rejec¬ 
tion of the Galaxy while minimising the number of modes 
lost from the analysis. Working with Pa ensures that we 
have projected out those modes of a which are (nearly) lo¬ 
calised in the region external to S 2 '. This procedure is sim¬ 
ilar to power spectrum estimation from the cut map, since 
projecting pixelised data into S 2 ' in real space has the effect 
of greatly amplifying the noise on the multipoles estimated 
from the cut map along those directions in multipole space 
which correspond to functions nearly localised in the cut. 
In practice, it is likely to be more efficient to represent the 
projected data vector Pa on a basis adapted to the region 
S 2 ', so that we work with the object U^a, where U is the 
(non-square) matrix whose columns are the eigenvectors of 
P with unit eigenvalue. The matrix U can easily be obtained 
from a singular value decomposition of P (Mortlock et al. 
2001). 


6 CONCLUSION 

Most current techniques for analysing CMB data on the 
sphere are based on the use of pixelised maps. The basis 
for a complementary approach, based on spherical harmonic 
coefficients of the intensity and polarization, has been de¬ 
rived here for the case of experiments that scan the full 
sky in circles. Our method offers a number of advantages, 
most notably the refined treatment of non-ideal beam ef¬ 
fects, the ability to handle correlated (low frequency) noise 
in an optimal, but efficient, manner, and the seamless way 
in which polarized data can be analysed alongside total in¬ 
tensity data. In principle, harmonic methods can be used to 
develop an end-to-end pipeline for the analysis of full-sky 
survey data, allowing the clean propagation of noise and 
other errors from the time-ordered data to the CMB power 
spectra. In practice, the methods described here are best 
regarded as being complementary to more standard map- 
based techniques rather than a replacement. Undoubtly, 
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there are analysis projects that are better suited to map- 
based techniques - typically those concerned with the sci¬ 
ence of local features in the maps, such as foregrounds. In 
addition, power spectrum estimation becomes cumbersome 
with purely harmonic methods if we have reason to question 
the foreground contamination of certain linear combinations 
of the spherical multipoles (e.g. those corresponding to fea¬ 
tures localised in the galactic plane). 

The biggest hurdle facing a practical implementation 
of the harmonic ‘map’-making method at the resolution de¬ 
manded by the upcoming satellite missions, is the need to 
solve the 0(1 ^ ax ) x O(Zmax) linear system in equation (28). 
Unlike conventional map-making in the presence of noise 
correlations, the problem lies not in the inversion of the noise 
covariance matrix, since we are working in a representation 
where this matrix is already very sparse. The inversion of 
A f l\l -1 A can be avoided by adopting iterative techniques 
with a block-diagonal preconditioner. The main computa¬ 
tional overhead arises from computing and storing the large, 
non-sparse matrix A and applying it to the sky multipoles 
and its transpose to the Fourier ring data. These operations 
can be performed very efficiently using fast Fourier trans¬ 
form techniques for the rather idealised case of constant lat¬ 
itude scanning, with rings uniformally spaced in azimuth 
and no variations in instrument properties during the mis¬ 
sion (Wandelt & Gorski 2001). Unfortunately, it does not 
appear that such techniques can be easily extended to more 
realistic scan strategies. An assessment of some of these nu¬ 
merical problems, together with a number of potential solu¬ 
tions, will be given in Mortlock et al. (in preparation). 
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APPENDIX A: MAIN BEAM VARIATION 
ACROSS A SPECTRAL BAND 

In this appendix we discuss the systematic error that arises 
from ignoring the variation of the beam multipoles with fre¬ 
quency across a spectral band when analysing microwave 
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data. This question does not appear to have received much 
attention so far in the CMB literature. For simplicity, we 
consider a single detector which is only sensitive to the total 
intensity. We take the beam profile to be an axisymmetric 
Gaussian, with frequency-dependent beam width cr(iz), so 
that for a(u) C 1 we have 

= SmO^^W,(u), (Al) 

where the scalar window function is 

Wi(v) = exp[— 1(1 + l)o 2 (v)/2}. (A2) 

Assuming that the beam is diffraction limited, we expect the 
frequency dependence of the beam width to go like a{v) = 
(vo/v)<jo, where <ro = u(t'o) is the beam width at the central 
frequency of the band. For l 1 this gives 

(^o), (A3) 

which also holds more generally for arbitrary beams which 
scale inversely with frequency along longitudes. 

If the CMB is the dominant physical component in the 
frequency band under consideration, the frequency spectrum 
of the brightness anisotropies in linear theory is 



500 1000 1500 2000 


l 

Figure Al. The systematic error in the recovered temperature 
power spectrum from neglecting the variation of the beam size 
across the spectral band for the 100 (solid lines), 143 (dashed 
lines), and 217 GHz (dash-dotted) Planck HFI channels. Also 
plotted is the random error in the C; from the combined instru¬ 
ment (white) noise from all detectors at the specified frequency, 
and from cosmic variance. We have assumed that all beams are 
axisymmetric with Gaussian profiles, and that the beam size is 
diffraction limited. 




dB[v, T ) 


dT 


(A4) 


for l > 0, where B(y, T) is the Planck function and T is the 
average CMB temperature over the sky. The which 

describe the dimensionless anisotropy in the thermodynamic 
temperature of the sky, are independent of frequency for the 
CMB. We further restrict attention to white noise, assume 
uniform coverage of the spin axis pointing over the entire sky, 
and ignore any systematic variation in focal plane geometry. 
Then, if we attempt to solve for the but ignore the 

variation of Wi(v) with frequency, the maximum-likelihood 
solution returns a biased estimate af lm ) which has the form 


-i 

a (ln 


— (1 + 5l)a\lm) + n (ln 


(AS) 


where the fractional bias 

/ Wi(y)(dB/dT)v{y)dv 
Wi(vo)f(dB/dT)v(v)dv 


(A6) 


and the dimensionless random errors have zero mean, 

and covariance 


(n(i )> = S W S mm' urVf 2 ( u 0 ). (A7) 

Here, w is the dimensionless weight per solid angle (Knox 
1995; see also Section 3). We can use the to estimate 
the CMB power spectrum Ci with the estimator 

Ci = 2T+1 Z l4m)| 2 - w-VfVo), (AS) 

| m | £ 

which would be unbiased if 5i were zero. The fractional sys¬ 
tematic error in our estimate of Ci is therefore 


{ACi/Ci) syst =6i(5i+2), (A9) 

while there is a fractional random error 

(A<7«/C,)rand = JzTfjl 1 + cr'w-'w,- 2 ^ „)]. (A10) 


The first term in brackets on the right-hand side of 
equation (A10) is cosmic variance. If instead we asked 
how well we could reconstruct the rotational-invariant 
71 |fij im )| 2 /(2? + 1) for our given realisation of the sky, the 
cosmic variance term would not be present. The latter sit¬ 
uation is the more relevant for mapping other astrophysical 
foreground components. 

In Fig. Al we plot the systematic and random errors 
on the CMB power spectrum for the 100, 143, and 217 GHz 
channels of the Planck HFI, using the predicted instrument 
specifications. The spectral filter v(v) is assumed to be a top 
hat with fractional width Av/vq = 0.33. It is clear that in 
the 100 GHz channel, the systematic error due to our ne¬ 
glect of the variation of beam width across the spectral band 
is non-negligible compared to the random error due to in¬ 
strument noise over a broad range of l. For an axisymmetric 
beam the effect can be easily accounted for by including the 
frequency-dependent window function Wi (n) in the integral 
over frequency, since the integral / af^ iy)b^ (lml) ( v)v{v ) dv 
would still factor into a part depending on the sky with 
indices l and m only, as required for subsequent analysis. 
However, for a non-axisymmetric beam where the variation 
of &d(im , )( iy ) cannot be reduced to a frequency-dependent 
factor with only an l index, and a constant part with l and 
m! indices, we can no longer integrate the variation of beam 
with frequency and still preserve the factorisation of the sky. 
Note that if there is only an effective beam asymmetry, aris¬ 
ing from skewing a symmetric beam with the temporal re¬ 
sponse of the instrument, the above comments do not apply 
since the frequency dependence of the effective beam will 
factor. 
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APPENDIX B: CONNECTING RING-SETS 
WITH MAPS 


In this appendix we establish the result given as equa¬ 
tion ( 66 ) in the main text, which relates the geometry of 
the ring-set to the observing time per solid angle on the sky. 

The first step is to replace the single summation 
Xq n |<min(! i') 011 the left-hand side of equation ( 66 ) with 
the double summation y\„. , , , Sinn y\ ,, . ,,,,,,, 

and to replace the labels l' by l" in the argument of the 
summation. Denoting the left-hand side of equation ( 66 ) by 

rll 1 


slmm' — ^ \ Svi"S m m' \/( 2 £ + l)(2l" + 1 ) 

l"n 

X d l mn (0)dn n (0)d l ns (a)d l ns{a)\, (Bl) 


where we have suppressed the limits on the summations. We 
now use the orthonormality of the spin-weighted harmonics, 


f —si ( I'm /) (^p) —si (Z // m) (^p) dflp — dmm / , (B2) 


to replace the Kronecker deltas in equation (Bl) with the 
integral of spin-weight — s harmonics over the sphere. The 
term (e p ) can be expressed in terms of D-matrices 

using equation (47): 

( e p) = Or, M, (B3) 

where we have also used d( nm ,(/I) = (-l) m ~ m ' d l m , m (/3) (e.g. 
Brink & Satcliler 1993), so that 

l"n 

x [(21” + 1 )d l rnn (0)d l ns (a)d l srn (6 p )\ 
x d l mn (e)d l rls (a)}dn p . (B4) 

If we now reverse the order of summation in equation (B4) 
using 

E E = E E • < B5 > 

Z // ^max(|m|,|s|) |n|^min(Z,Z // ) |n|^Z Z // ^max(|ra|,|n|,|s|) 

the summation over l" can be performed with the Ponzano- 
Regge sum rule (e.g. Varshalovich, Moskalev & Khersonskii 
1988; equation 21, p. 89): 

^ (2/" + l)d l mn (9)d l ns (a)d l srn (9 p ) 

l ,, ’^max(\m\ ,|n|, |s|) 

= — cosfmdi + nd 2 + sd 3 ), (B6) 

w y/B v 

where B p is given by equation (61), and the angles <5i, S 2 , 
and S 3 are defined by the SO (3) composition 

D(Si, 9, 0)D(52, a, 83 ) = D( 0, — 9 P , 0). (B7) 

(The angles are given explicitly by Varshalovich et al. 1988, 
but it is straightforward to show that their definitions are 
equivalent to the implicit definition given here.) Finally, we 
can perform the remaining sum over n in equation (B4) by 
writing cos(m<5i +n<52 + sd 3 ) as the real part of exp[— i(m8i + 


71S2 + 5 ^ 3 )], and combining the complex exponentials with 
the remaining d-functions to get 

y d l mn ( 6 )d l ns (a) cos(m<5i + n<5 2 + sds) 

|n|^Z 

= SR ^ D l rnn ( 5 i, 9 , 0)D l ns (82, a, S 3 ) 

= d l sm (0 P7 ), (B 8 ) 

where we have used equation (B7) in the last equality. The 
remaining terms in equation (B4) combine with d l sm ( 9 p , ) to 
give -sl''(; m )(e p ) on using equation (B3). Our final result is 

2 -^d. s Sr im) (e p ) - s Y (l , m/) (e p ) dBp, (B9) 

n V b p 

which establishes equation ( 66 ). 

This paper has been produced using the Royal Astronomical 
Society/Blackwell Science style file. 
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